Scale‐dependent environmental effects on phenotypic distributions in Heliconius butterflies

Abstract Identifying the relative importance of different mechanisms responsible for the emergence and maintenance of phenotypic diversity can be challenging, as multiple selective pressures and stochastic events are involved in these processes. Therefore, testing how environmental conditions shape the distribution of phenotypes can offer important insights on local adaptation, divergence, and speciation. The red‐yellow Müllerian mimicry ring of Heliconius butterflies exhibits a wide diversity of color patterns across the Neotropics and is involved in multiple hybrid zones, making it a powerful system to investigate environmental drivers of phenotypic distributions. Using the distantly related Heliconius erato and Heliconius melpomene co‐mimics and a multiscale distribution approach, we investigated whether distinct phenotypes of these species are associated with different environmental conditions. We show that Heliconius red‐yellow phenotypic distribution is strongly driven by environmental gradients (especially thermal and precipitation variables), but that phenotype and environment associations vary with spatial scale. While co‐mimics are usually predicted to occur in similar environments at large spatial scales, patterns at local scales are not always consistent (i.e., different variables are best predictors of phenotypic occurrence in different locations) or congruent (i.e., co‐mimics show distinct associations with environment). We suggest that large‐scale analyses are important for identifying how environmental factors shape broad mimetic phenotypic distributions, but that local studies are essential to understand the context‐dependent biotic, abiotic, and historical mechanisms driving finer‐scale phenotypic transitions.


| INTRODUC TI ON
Investigating how environmental conditions influence the distribution of species and their populations can provide important insights into the generation and maintenance of phenotypic diversity (Edelaar, 2018;Gould & Johnston, 1972;Schluter, 2001). Indeed, numerous studies have documented local phenotypic adaptation in response to habitat heterogeneity (e.g., Barrett et al., 2019;Boncoraglio & Saino, 2007;Hegna et al., 2013). However, it has long been recognized that a combination of selective forces and stochastic events may influence the spatial distribution of phenotypes (Amézquita et al., 2009;McLean & Stuart-Fox, 2014). Accordingly, disentangling the factors promoting phenotypic diversity remains a complex task, as these mechanisms can directly or indirectly influence each other (Gray & McKinnon, 2007;Nosil, 2012).
Heliconius butterflies represent an excellent system for exploring phenotypic diversity and identifying the environmental factors influencing phenotypic distributions. Distantly related Heliconius species (along with other butterfly and moth genera) have converged on similar aposematic wing colors, forming a handful of local mimicry rings (Birskis- Barros et al., 2021;Mallet, 1999;Papageorgis, 1975;Turner, 1971). Interestingly, the presence of various aposematic rings in a single area contradicts predictions from the Müllerian mimicry theory, which posits that local phenotypic convergence resulting in shared warning patterns should reduce the cost of predator learning (Joron & Mallet, 1998;Mallet, 1993;Sherratt, 2008). This system also presents a remarkable variation of wing color patterns within a single mimicry ring, such that populations of a species characterized by specific colors (such as red, yellow, and black) can have very different patterns, resulting in patchworks of distinct phenotypes across the range of a species (Brown Jr. et al., 1974;Papageorgis, 1975;Turner, 1971).
While previous research has shown that wing color in Heliconius populations are under strong natural and sexual selection (Dell'Aglio et al., 2016, 2018Jiggins et al., 2001Jiggins et al., , 2004Mallet & Barton, 1989;Merrill et al., 2012Merrill et al., , 2014, the environmental factors driving the spatial distribution of color phenotypes among and within mimicry rings remain unclear (Gompert et al., 2011;Joron, 2005;Joron & Mallet, 1998;Mallet, 1993). More specifically, although some studies have tested the influence of environmental gradients on Heliconius phenotypes (Arias et al., 2008;Blum, 2008;Jiggins et al., 1996;Rosser et al., 2014;Thurman et al., 2019), most large-scale distribution maps are based solely on occurrence points, without environmental variables being used as predictors (Brown Jr. et al., 1974;Brown Jr., 1982;Rosser et al., 2012;Turner, 1971). Additionally, previous distribution modeling of Heliconius species has not accounted for intraspecific phenotypic variation (Rueda-M et al., 2021). This is relevant because distinct phenotypes usually have varying fitness advantage depending on the location (Arias et al., 2008;Blum, 2008), leading to differences in their geographical distributions. Finally, no study has investigated the associations among environmental variables and Heliconius phenotypes across distinct spatial scales, despite a growing recognition that the relative influence of environmental predictors may be scale-dependent (Bunnell & Huggard, 1999;Sandel & Smith, 2009;Wiens & Bachelet, 2010).
In this study, we estimated the spatial distribution of Heliconius phenotypes within a single mimicry ring, using two distantly related and non-hybridizing co-mimics, Heliconius erato and Heliconius melpomene. These two species are part of the red-yellow mimicry ring [red and yellow group according to Turner, 1971], which shows a remarkable phenotypic variation in red, yellow and orange wing elements across the Neotropics (Figure 1a). We used a multiscale approach including a large-scale analysis covering the entire distribution of the co-mimics, as well as a local-scale analysis focusing on hybrid zones within the Brazilian Amazon Forest. The Amazon region contains some of the most biodiverse terrestrial ecosystems (Da Silva et al., 2005;Myers et al., 2000), which may be partially attributed to high variation in precipitation, seasonality, and temperature (Hijmans et al., 2005). Our local-scale sampling covers regions explored by the naturalist Henry Walter Bates over 150 years ago.
Notably, Bates postulated that the geographical turnover in H. melpomene phenotypes he observed during his explorations was the result of variation in humidity, vegetation, and soil types (Bates, 1862).
By analyzing the effect of different environmental variables on the distribution of the distinct phenotypes belonging to the redyellow mimicry ring ( Figure 1a) at multiple scales, we aim to answer: (i) Is there a general relationship between Heliconius wing color phenotypes and environmental variables across the Neotropics?
(ii) How congruent are the predictive variables between co-mimics?
(iii) If phenotype-environment associations exist, are they consistent across distinct spatial scales? In addition, we aim to find support for two hypotheses based on different mechanisms of selection: (1) The Signaling hypothesis, which predicts that phenotypes with larger melanized wing areas are associated with open-canopy forests and higher solar irradiation, since darker aposematic patterns are more conspicuous against lighter environment backgrounds; accordingly, we expect that paler phenotypes should be found in denser forests, as these patterns better advertise aposematism in shaded environments (Arenas et al., 2014;Osorio & Vorobyev, 2005;Rojas et al., 2018;Stevens, 2007); and (2) The Thermoregulation hypothesis, which predicts that darker phenotypes are associated with cooler environments with lower solar irradiation, as darker wings provide greater efficiency for thermoregulation; accordingly, we expect that less melanized phenotypes should be found in warmer habitats with higher solar irradiation, as these would be more efficient at avoiding overheating (Clusella-Trullas et al., 2007;Clusella-Trullas & Nielsen, 2020;Hegna et al., 2013;Van Dyck & Matthysen, 1998).

| Large-scale analysis (distribution models)
We used Species Distribution Modeling (SDM; Guisan & Thuiller, 2005) to predict the distribution of H. erato and H. melpomene phenotypes within the red-yellow Müllerian mimicry ring (Figure 1a).
SDM is an approach that combines known coordinate occurrences with information about the environmental tolerances of a species to map suitable habitats (Austin, 2002;Elith & Leathwick, 2009).
We sampled butterflies using entomological nets from December 2016 to June 2017 and from July to October 2018 at ~25 km intervals across four transects in the Brazilian Amazon, totaling 77 sites distributed across ~1905 km (Figure 1b; Table S1). Transects were selected to capture a wide array of ecosystems ( Figure S1) and variation in environmental and physical features across the Amazon, as well as to include regions containing hybrid zones within the two studied species. We also included samples from the collection of the

Laboratory of Ecology and Systematics of Pollinators and Predators
(LESPP-Federal University of Maranhão) and data used in Rosser et al. (2012), available at https://helic onius -maps.github.io/. For each species, subspecies belonging to the same phenotype were considered as a single group to produce a predicted distribution map (Table S2). In total, we used 3403 occurrence points after filtering for duplicated records, uncertainty or incomplete taxon identification, and grids without environmental data ( Figure S2).
Our environmental dataset comprised 31 variables selected based on their potential direct and indirect effects on Heliconius ecology and evolution (Table S3). Selected variables were prepared using QGIS (QGIS Association, 2020) with a 1 km resolution. To avoid collinearity and model overfitting, we performed Pearson's correlation and selected variables using a cut-off of r < .7 (Dormann et al., 2013). Our final dataset comprised 12 environmental variables, including climate, soil, vegetation, wind, and solar parameters (Table S3).
We performed distribution modeling using the Maximum Entropy Method (MaxEnt; Phillips et al., 2006) with the SDMtune package (Vignali et al., 2020) in R (R Core Team, 2019). Maxent is a machine learning method that outperforms other distributional modeling techniques based on different metrics that assess the agreement between records and predictions (Elith et al., 2006;Phillips et al., 2006;Wisz et al., 2008). Maxent does not require true absence data, which is rarely available (Elith et al., 2006;Phillips et al., 2006), and its results can be interpreted as habitat suitability for taxa based on a set of environmental predictors (Phillips et al., 2006). For background data, we randomly created 10,000 points covering the entire distribution of H. erato and H. melpomene (Barbet-Massin et al., 2012).
We randomly selected 75% of presence and background data for training, and 25% for testing models' predictability. Additionally, we used a 10-fold cross-validation test as a validation strategy (Braunisch & Suchant, 2010) and performed the fine-tuning strategy to select Maxent's best parameters (Vignali et al., 2020; Table S4).
This method allowed us to achieve the best predictive power, while avoiding model overfitting (Hallgren et al., 2019;Radosavljevic & Anderson, 2014).  metrics, considering an AUC > 0.7 and a TSS > 0.5 as a non-random distribution model with respect to the predictors (Allouche et al., 2006;Elith et al., 2006). To evaluate the importance of each environmental variable as a model's predictor, we used the permutation importance metric. This was calculated by permuting one variable at a time (using 10 repetitions) and computing the decrease in the model's AUC metric. Excluded variables that mostly decreased AUC values were considered more important to the model (Vignali et al., 2020).
We converted predicted distributions into binary maps using the maximum sensitivity plus specificity as threshold (Liu et al., 2013(Liu et al., , 2016. The raster R package (Hijmans & van Etten, 2020) was used to calculate the intersection area in the distribution of co-mimics and Logistic Maxent response curves were used to evaluate the probability of a phenotype being present in a locality (Merow et al., 2013).
Finally, we overlapped the predicted distribution of each phenotype with a map of the world's biomes (Dinerstein et al., 2017) to evaluate whether changes in phenotypes are associated with biome boundaries. As biomes represent global-scale ecosystems characterized by unique sets of environmental conditions (Hoekstra et al., 2004;Olson et al., 2001), we expected to find a pattern of changes in phenotypes across biome boundaries.

| Local-scale analysis
To test the effects of local-scale environmental gradients on phe-  Table S5) and performed a Principal Component Analysis (PCA) with the patternize R package . To assess differences in black patterning among phenotypes, we estimated the wing surface area containing black scales in the different clusters using one-way ANOVA and Tukey's tests ( Figure S3).
For our local-scale analysis, we focused on hybrid zones involving pure phenotypes that presented distinct black patterning (higher melanin x lower melanin). Using this approach, we were able to test opposing hypotheses about mechanisms of selection from visual signaling versus thermoregulation. Importantly, we could not perform this analysis in our large-scale approach because we did not have a scenario in which there were clear transitions of morphs showing a significant difference in black patterning ( Figure S3).
We performed cline analyses using a Bayesian approach implemented in the bahz R package (Thurman, 2019). Geographic coordinates were transformed using linear regression to obtain onedimensional transects, adapting a method described in Thurman et al. (2019). The "new" coordinates were set as the closest point within the linear transect to the true location (Figure 1e), and distances were calculated using two approaches: (1) relative to the westernmost or northmost sampling site for each species, which were not necessarily the same as we were not able to collect both species in all the sites and (2) (Mallet, 1986) and/or sparse yellow scales on the dorsal hindwing.
We used cline results to overlay hybrid zones on a map of terrestrial ecoregions (Dinerstein et al., 2017) to test whether changes in phenotypes are associated with ecoregion boundaries. Ecoregions are defined as regional-scale ecosystems nested within biomes (Hoekstra et al., 2004;Olson et al., 2001); as such, ecoregions represent ecological boundaries at finer scales when compared with biomes.
Finally, to directly evaluate whether changes in phenotypic frequencies at local scales are associated with environmental variables, we performed spearman correlations using seven environmental parameters, including climate, solar and vegetation variables (Table S3). Only variables with pairwise correlations lower than 0.7 (to avoid collinearity; Dormann et al., 2013) and that were found to be important in the best-fitting distribution models in our large-scale analysis were selected.

| Abiotic variables predict large-scale phenotypic distributions
We found that environmental variables play an important role in driving the large-scale distribution of phenotypes within the Heliconius red-yellow mimicry ring, as indicated by the high values of our models' predictive performance (AUC > 0.8 and TSS > 0.5, Table S6). The The exception was the dennis phenotype, in which annual mean temperature was the most important predictor for H. erato (29.04%), and precipitation of the warmest quarter (33.57%) for H. melpomene.
However, temperature seasonality is still among the most important predictors for dennis models, being the second highest contribution for H. erato (28.23%) and the fourth highest contribution for H. melpomene (8.76%; Figure 2).
The rank order of contributions by environmental predictors showed that thermal and precipitation seasonality are particularly important for predicting phenotypic distributions of H. erato and H. melpomene co-mimics, although precipitation seasonality had lower relative importance for PHYB. Percentage of tree cover and solar irradiation (GHI) were also important for predicting distributions, although less important for the dennis phenotype ( Figure 2). In general, soil parameters and wind speed had low contributions in most models.
Response curves (i.e., the probability of a phenotypic presence as a function of an environmental gradient) indicate that environmental predictors are more similar between co-mimics than among phenotypes within the same species (Figures 3 and S4). For instance, two thermal variables (temperature seasonality and annual mean temperature) showed distinct responses among phenotypes but generally similar optima values and response curves between co-mimics in three out of the four co-mimetic phenotypes. The exception was the PHYB phenotype, in which the phenotypic optimum between co-mimics were highly different. Our maps also showed considerable overlap (>47%) in the predicted distributions of co-mimics, except for the PHYB phenotype (22.6%; Figure 4). Interestingly, even though several environmental variables contributed to explaining the large-scale distribution of Heliconius phenotypes, we did not find abrupt changes in phenotype composition across biome boundaries.
In fact, our models suggest that all phenotypes occur in more than one biome (Figure 4). Notably, the postman phenotype was associated with warmer environments and lower temperature seasonality when compared with the dennis and dennis-ray phenotypes ( Figures 3 and S4).

| Effects of abiotic gradients on local phenotypic transitions
Consistent with the broad scale patterns, we found that local phenotypic transitions are not associated with ecoregion boundaries.
Nonetheless, in contrast to our large-scale results, our local-scale analyses revealed incongruencies in co-mimic relationships with environmental gradients and inconsistent environmental effects in different regions of the Amazon.
In transition 1 (Figure 1c1), we observed a similar shape and location for hybrid zones of both H. erato and H. melpomene. This transition followed the change from dennis-ray (a lower melanin morph- Figure S3) to postman (a higher melanin morph- Figure S3).
The predicted center of the H. erato hybrid zone is 310.78 km from the westernmost collection site and has a predicted width of 173.68 km (Figure 5a; Table S7). H. melpomene has a steeper hybrid zone with a center at 355.68 km from the westernmost collection site and has a width of 93.74 km (Figure 5a; Table S7). Taking into consideration the common westernmost site for both co-mimics, the center of the H. melpomene hybrid zone is at 331.07 km and has a width of 92.66 km (Figures 5a and S5; Table S7). Both hybrid zones are entirely located within a single ecoregion (Figure 5a), indicating that changes in ecoregions are not associated with this phenotypic transition. However, changes in specific environmental variables show congruent associations with phenotypic transitions between co-mimics. We found that the frequency change in the postman phenotype across the hybrid zone for both species and environmental variables indicate that compared with the lower melanin dennisray phenotype, the higher melanin postman phenotype is associated with higher mean temperatures, higher precipitation seasonality, and higher solar irradiation, as well as lower annual mean precipitation, lower temperature seasonality, and lower vegetation density ( Figure 6a; Table S8).
Transition 2 (Figure 1c2) occurs in the same region as transition 1 and follows the change from postman phenotype (a higher melanin morph- Figure S3) to the PHYB phenotype (a lower melanin morph- Figure S3). Within this phenotypic transition, there is incongruence between the co-mimic hybrid zones. The predicted center of the H. erato hybrid zone is at 208.10 km from the westernmost site and has a predicted width of 194.08 km, crossing three different ecoregions (Figure 5b; Table S7) Table S8). However, within these phenotypic transitions we found that, compared with the lower melanin PHYB phenotype, the higher melanin postman phenotype is associated with lower annual mean temperature, temperature seasonality, and precipitation seasonality ( Figure 6b; Table S8).

| DISCUSS ION
Our multiscale approach revealed that phenotypic distributions were driven by environmental variables (predominantly climatic variables associated to heat and pluviometry seasonality), with congruent patterns between most H. erato and H. melpomene co-mimics across the Neotropics. However, we found more nuanced patterns at finer spatial scales, with phenotype-environment associations being congruent between co-mimics in some hybrid zones and in-

| Phenotypic distribution across the neotropics
At a broad spatial scale, phenotypic distributions were mostly explained by temperature seasonality and to a lesser extent by precipitation seasonality. It is likely that seasonality directly affects  Freq. postman is currently no evidence linking specific color patterns within H.
erato and H. melpomene to higher or lower tolerance of heat and/ or dry conditions. Seasonality may also have an indirect effect on Heliconius distributions due to effects on hostplant distributions and abundance, as well as on Heliconius hostplant preference range (Benson, 1978;Jorge et al., 2011;Merrill et al., 2013;Rodrigues & Moreira, 2002. Even if the abiotic conditions are suitable for a given Heliconius phenotype, populations may not persist locally because of a lack of suitable hostplants (Araújo & Luoto, 2007;Lemoine, 2015).
Our large-scale findings also indicate that the postman phenotype is predicted to occur in warmer regions with more constant temperatures compared with the dennis-ray and dennis phenotypes.
Interestingly, we found no congruent patterns for the postman with a yellow hindwing bar phenotype (PHYB), with distinct distributions and associations with environmental gradients between co-mimics. PHYB is more restricted to the latter. Importantly, differences in large-scale phenotypic distributions could also be a consequence of distinct colonization and/or extinction histories (Brower, 1996;Joron & Mallet, 1998;Mallet et al., 1996;Quek et al., 2010;Turner & Mallet, 1996).

| Local phenotypic distributions
Environmental variables were also important in predicting local distributional patterns in most of the transition zones we investigated. However, at the local level there were differences in the congruence of environmental effects between co-mimics and in the consistency of phenotype-environment associations across distinct hybrid zones. In the two hybrid zones located in the Eastern Amazon, we found that co-mimics share relatively similar shaped transitions (i.e., curves showing phenotypic frequency change along the transect) and that those transitions were associated with similar environmental variables. However, the specific environmental features differed among these hybrid zones. The dennis-ray to postman transition shows greater congruence between co-mimics than the postman to PHYB transition. In contrast, there was no congruence between comimics in the third hybrid zone located in the center of the Amazon. This is true for both the phenotypic frequency curve shapes and the associations with different environmental variables.
A number of factors may be interacting and are likely important in shaping local-scale distributions of species and populations. These include local selective pressures (Lenormand et al., 2009;McLean & Stuart-Fox, 2014), geological history (Hall & Harvey, 2002), shared ancestry (Losos, 2008;Wiens & Graham, 2005), different diversification histories (Quek et al., 2010), and stochastic events causing genetic drift (Mallet, 1993(Mallet, , 2010Sherratt, 2006). These factors are probably interacting across the area of the Amazon basin that we studied. At the population level, frequency-dependent selection (FDS) by predators, a selection regime in which the fitness of a local phenotype increases with its frequency (reviewed in Ruxton et al., 2004;Sherratt, 2008), is thought to be among the most important forces influencing aposematic pattern distributions in Heliconius (Finkbeiner et al., 2014;Kapan, 2001;Mallet & Barton, 1989). The magnitude of FDS is influenced by the presence of different predators (Nokelainen et al., 2014;Willmott et al., 2017), their distinct visual systems (Dell'Aglio et al., 2018), and innate differences in their learning and memory abilities (Endler & Mappes, 2004;Mappes et al., 2005;Speed & Turner, 1999), as well as the prey density (Kapan, 2001) and the dispersal of the butterflies (Barton & Hewitt, 1985. In this light, one might predict that selection on wing patterns that are more effective visual signals in a given environment might explain local phenotypic transitions. However, of the three transitions zones that we assayed, only one is consistent with this hypothesis. The postman phenotype possesses higher percentage of wings covered with black scales (proxy for the amount of melanin) than the dennis-ray phenotype and would be a better aposematic signal in warmer, brighter and open habitats (Arenas et al., 2014;Osorio & Vorobyev, 2005;Rojas et al., 2018;Stevens, 2007). The transition from dennis-ray to postman phenotype, located in Eastern Amazon, is consistent with this expectation.
Indeed, molecular variation around genomic regions responsible of red color patterns suggests that the dennis-rayed phenotype arose recently in both H. erato and H. melpomene and spread across the Amazon basin replacing the more ancestral postman phenotype (Hines et al., 2011). The hypothesized spread of this phenotype could be due to its better signaling ability in moist Amazonian forests. This phenotype is part of a much larger mimicry ring composed of over a dozen different butterfly and moth species (Papageorgis, 1975;Turner, 1971). Thus, a combination of signal efficiency and an experienced predator community could be responsible for this association. The genetic evidence for widespread phenotypic replacement remains speculative and our phenotypic data is only correlational.
Thus, we cannot exclude the action of other selective pressures or random processes (Cuthill et al., 2017;Postema et al., 2022;Schluter et al., 1991) for the contemporary phenotypic distributional patterns. Moreover, in a hybrid zone approximately 1100 km away and involving the same phenotypes, there was no compelling support for the signaling hypothesis. Here, the H. erato and H. melpomene transition zones were disjunct. The H. erato zone was centered on the Amazon River, which is a potent, albeit often semipermeable, dispersal barrier for many species (Godinho & da Silva, 2018;Hayes & Sewlal, 2004;Lynch Alfaro et al., 2015;Rosser et al., 2020;Smith et al., 2014). In contrast, the transition between H. melpomene phenotypes was located approximately 120 km south of the Amazon River and largely uncorrelated with environmental variables. These discordances probably reflect different dispersal histories (Quek et al., 2010;Van Belleghem, Baquero, et al., 2018).
Given the broadscale association of phenotypes with environmental variables, it also seems likely that differences in the ability for different wing pattern phenotypes to thermoregulate may help explain local distributional patterns. However, again the results of our local scale analysis were mixed, with only one transition zone being roughly consistent with differences in thermoregulation.
Here, the PHYB phenotype (i.e., postman with a hindwing yellow bar) is mostly associated with warmer environments, where its paler wings (compared with the postman phenotype) might provide a more efficient way to avoid overheating (Clusella-Trullas et al., 2007;Clusella-Trullas & Nielsen, 2020;Hegna et al., 2013;Van Dyck & Matthysen, 1998 A consistent theme in prior hybrid zone work is the suggestion of a role of the environment in shaping distributional patterns (Arias et al., 2008;Benson, 1982;Jiggins et al., 1996;Muñoz et al., 2010).
Accordingly, even if other factors show to have greater effects on hybrid zone locations (such as the balance of predation and dispersal), direct and/or indirect habitat effects cannot be completely dismissed (Mallet, 1993;Mallet & Barton, 1989). For example, Blum (2008), showed that color pattern genotype and phenotype frequencies corresponded to land cover differences in a French Guiana H. erato hybrid zone, supporting the hypothesis of visual signaling efficiency for mate choice. On the other hand, Mallet and Barton (1989) suggested that patterns of hybrid zones of H. erato and H. melpomene in Peru are primarily explained by predation; however, they did not discard subspecies adaptations to local ecological conditions. Thus, although the specific environmental mechanisms and the strength of their effects vary across localities, there is abundant evidence for the environment acting as an important driver of local patterns of phenotypic distribution.
Heliconius phenotypic variation associated with habitat characteristics can also be explored under the context of spatial-temporal dynamics, as hybrid zones can move in response to environmental gradients (Barton, 1979;Barton & Hewitt, 1985;Buggs, 2007;Mallet, 1986Mallet, , 1993Wielstra, 2019). Blum (2002), for example, suggested that deforestation could be driving the movement of a H. erato hybrid zone in Panama, which was characterized 17 years earlier (Mallet, 1986). This zone between a postman (H. e. hydara) and a PHYB (H. e. demophoon) subspecies continued to move over the next 15 years, but at a slower rate . Thurman et al. did not find a strong association between movement and deforestation, suggesting that the genetically dominant postman is replacing the PHYB phenotype through a phenomenon referred to as "dominance-drive" (see Mallet, 1986 in the region of very high precipitation. Therefore, the majority of emerging data, including this study, supports the role of the environment in shaping patterns of underlying genetic and associated phenotypic variation, even though other factors can interact and/ or mask its effect. One of the most compelling aspects of studying the distribution of phenotypes are the implications for local adaptation and ecological speciation (Nosil, 2012).   formal analysis (equal); investigation (equal); project administration (equal); resources (lead); supervision (lead); writing -original draft (equal); writing -review and editing (equal). Jonathan S.

ACK N OWLED G M ENTS
This study was financed by the Coordenação de Aperfeiçoamento

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.

FU N D I N G I N FO R M ATI O N
This study was financed by the Coordenação de Aperfeiçoamento